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ABSTRACT 


It  is  the  purpose  of  this  paper  to  develop  a  useful  mathematical 
model  of  ASW  aircraft  availability.      The  increasing  emphasis  of  systems 
studies  dictates  the  use  of  accurate  and  representative  models  of  the  ASW 
systems.     At  present,    many  studies  are  using  essentially  the  same  models 
developed  during  World  War  II.      This  paper  is  an  attempt  to  make  use  of 
advanced  theory  in  a  more  powerful  and  flexible  model  and  to  make  the  use 
of  the  model  practical  and  verifiable. 

The  writer  adapted  the  time  homogeneous  bivariate  model  as  de- 
veloped by  F.    C.    Collins.      This  is  a  discrete  time  Markov  process  with 
a  stochastic  matrix  of  transition  probabilities  wherein  the  maintenance 
process  is  modeled  as  a  pulsed  input  multiple  server  queue. 

The  model  was  programmed  in  FORTRAN  63  on  the  CDC  l604and 
then  modified  to  allow  for  variability,  in  the  input  parameters.      Other 
modifications  include  an  increase  in  the  size  of  the  model  to  accommodate 
a  16-air craft  squadron,    the  largest  ASW  squadron  at  present,    and  an 
explicit  form  solution  to  the  maintenance  queueing  equations. 
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requiring  maintenance  when  recovered  by  the 
carrier 

II      (m)  the  probability  that  of  a   a/c  flying   m    will  enter 
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1.       INTRODUCTION 

The  threat  to  freedom  of  the  seas  posed  by  the  vast  Soviet  sub- 
marine fleet  is  perhaps  the  most  thorny  problem  facing  the  U.  S.    Navy 
today.      Two  world  wars  have  produced  Pyrrhic  victories  over  limited 
submarine  fleets.      During  the  Second  World  War  operations  analysis 
was  born  into  the  Navy  to  aid  in  the  defeat  of  the  German  submarine. 
The  classic  antisubmarine  warfare  (ASW)  analyses  and  models  developed 
by  Morse  [2]   and  Koopmans  [3]   are  still  being  used  today,    over  two 
decades  later,    in  most  of  the  ASW  study  efforts  for  the  Navy. 

These  early  ASW  analyses  assumed  a  given  level  of  search  effort 
available  and  directly  evaluated  the  probability  that  an  ASW  subsystem 
could  detect  and/or  kill  a  submarine.      This  assumption  is  not  only 
logical  to  make  the  problem  tractable,    but  also  practical  since  no 
immediate  changes  in  ASW  force  levels  could  be  expected.     Moreover, 
the  studies  were  conducted  during  the  war,    not  before  it  started.      It  is 
the  purpose  of  this  paper  to  present  a  probabilistic  model  to  describe 
the  available  effort.      Such  a  model  can  be  used  to  sharpen  the  estimates 
of  the  effectiveness  of  an  ASW  subsystem  and  to  study  the  charac- 
teristics of  the  associated  support  system. 

Naturally,    the  current  study  plays  an  important  but  limited  role  in 
the  overall  problem  of  designing  an  entire  ASW  system.      The  diffi- 
culties involved  in  such  a  specification  are  legion.      First  and  foremost 
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is  the  quantification  of  the  ASW  mission  in  denying  the  enemy  the 
effective  use  of  his  submarines.      Currently,    the  probability  of  detecting 
and/or  killing  submarines  is  used  as  the  measure  of  effectiveness  of 
the  mission,    and  it  appears  that  a  more  encompassing  one  has  not  been 
developed.      Second,    the  specification  of  an  ASW  force  level  to  counter 
a  given  threat  has  many  inherent  subjective  elements.      These  are  due  to 
the  existing  historical  bias  in  predicting  the  conduct  of  a  future  ASW  war 
with  an  enemy,    particularly  one  who  has  never  before  used  a  large 
submarine  force  in  its  military  operations.      The  reader  can  imagine 
why  merely  defining  terms  such  as  "threat"  and  "effective  counter" 
becomes  quite  difficult. 

Thus,    there  is  a  need  to  investigate  the  levels  of  search  effort 
specified.      This  may  require  acceptable  models  to  measure  the  availa- 
bility of  effort,    its  effectiveness,    and  determine  the  logistic  support 
required  for  any  level  of  available  effort.      Specifically,    the  ASW  sub- 
system to  be  modeled  is  the  carrier -based  aircraft,    although  the  model 
is  adaptable  to  other  systems. 

The  method  of  investigating  the  demand  for  ASW  carrier  a/c  will 
assume  that  the  desired  number  of  a/c  on  station  is  known  as  an  input 
parameter.      The  support  required  to  achieve  this  measure  of  available 
effort  depends  upon  maintenance  space,    manpower,    and  supply. 
Generally,    we  shall  consider  how  an  ASW  carrier  supports  this  number 
of  a/c  on  station  with  the  present  or  proposed  number  of  a/c  embarked 
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on  the  carrier.      The  parametric  input  can  be  subjected  to  sensitivity- 
analyses. 

The  operational  commander  of  the  ASW  force  launches  the  desired 
number  of  a/c  on  station  to  screen,    search,    or  actively  prosecute  a 
submarine  contact.      Each  a/c  is  relieved  on  station.      Each  such  relief 
requires  the  launching  of  another  a/c  prior  to  the  recovery  of  the  initial 
a/c.      The  returning  a/c  must  receive  varying  degrees  of  maintenance 
and  requires  refueling  and  rearming.      This  cycle  continues  until  the 
mission  is  completed.      Loss  of  a/c  due  to  accidents,    insufficient  supply, 
and  lack  of  repair  capability  cause  deviations  in  this  procedure.      Naval 
operations  involve  the  interaction  of  many  quantities  which  are  random 
in  nature.     Not  all  can  be  considered  in  a  tractable  mathematical  model. 
Some  quantities  which  are  important  are  omitted.      One  example  is  the 
length  of  each  cycle  time,    which  is  assumed  to  be  a  constant  value. 
Including  variables  of  this  nature  incurs  unnecessary  mathematical 
complication.      It  is  hoped  that  adequacy  of  the  model  can  be  measured 
by  using  fleet  data  available  from  the  Fleet  ASW  Data  Analysis  Program 
(FADAP). 

Collins  [5]   describes  a  bivariate  Markov  model  for  airborne  early 
warning  (AEW)  and  combat  air  patrol  (CAP)  jet  a/c  operating  in  an 
attack  carrier  force.      This  model  is  used  to  evaluate  the  probability  of 
maintaining  a  fixed  requirement  of  a/c  on  station  as  a  measure  of  ef- 
fectiveness of  the  system.      It  has  subsequently  been  used  in  a  larger 
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attack  force  study  for  the  Navy.      The  model  computes  the  probabilities 
of  the  number  of  a/c  on  station  and  in  or  awaiting  maintenance  at  any 
given  launch  period.      The  comparable  ASW  problem  differs  in  the  following 
aspects: 

1.  Type,    range,    and  speed  of  a/c; 

2.  The  variable  number  of  a/c  required  for  mission; 

3.  Attrition  due  to  accidents  and  supply  failures; 

4.  The  greater  number  of  ASW  a/c. 

It  was  decided  to  use  the  Collins'  model  with  appropriate  modification. 
For  immediate  reference,    the  mathematical  content  of  the  model  will 
be  repeated  herein. 

In  order  to  incorporate  these  modifications,    it  was  necessary  to 
spend  some  time  reprogramming  on  the  CDC  1604  digital  computer  in 
FORTRAN  63,    the  CDC  version  of  the  IBM  FORTRAN  IV.      The  original 
program  was  not  readily  available  and  was  written  in  an  early  assembler 
language.      Moreover,    the  numerical  analysis  was  not  sufficiently  sharp 
to  handle  the  larger  input  values.      Also,    double  precision  (two  computer 
words  instead  of  one)  arithmetic  was  required  in  one  subroutine  for  an 
accurate  explicit  solution  to  the  maintenance  queueing  equations  (see 
Appendix  I).      This  effected  a  50%  decrease  in  the  computer  time 
required  for  developing  a  matrix  of  transition  probabilities. 

Following  this  introduction,    section  2  contains  a  brief  description  of 
the  operational  problems  involved  and  the  assumptions  made.     A  brief 
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description  of  Markov  chains  and  the  mathematical  model  are  presented 
in  section  3.      The  details  for  computing  the  matrix  of  transition  proba- 
bilities are  given  in  section  4.      General  employment  of  the  model 
follows.      The  appendices  include  the  solution  mentioned  on  the  pre- 
ceding page,    a  logical  flow  diagram  of  the  program,    a  copy  of  the 
program,    and  some  sample  results. 
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2.       ASSUMPTIONS 

The  real-world  employment  of  carrier  a/c  is  cyclic  in  nature,    and 
the  present  state  of  any  given  a/c  (i.  e.  ,    flying,    in  or  awaiting  mainte- 
nance) depends  largely  on  what  the  previous  state  was.      This  fact 
suggests  that  a  Markovian  assumption  can  logically  be  made  for  the 
a/c  transition  probabilities.      In  the  search  phase,    a/c  may  or  may  not 
relieve  on  station;  but,    in  any  part  of  the  contact  investigation  phase, 
relief  on  station  will  be  made.      To  insure  full  screening  and  mission 
coverage,    a/c  will  relieve  on  station. 

The  question  of  resupply  during  an  operation  depends  primarily  on 
the  availability  of  carrier  on-board  delivery  (COD).      This  depends  on 
the  geographical  location  and  the  mission  (convoy  protection,    strike- 
force  protection,    hunter -killer  operation,    etc.).     In  practice,    resupply 
is  not  anticipated  within  a  week's  period,    and  around-the-clock  oper- 
ations have  continued  for  two  weeks  without  resupply. 

Standard  maintenance  procedures  aboard  carriers  preclude  major 
maintenance  on  the  flight  deck.      It  will  be  assumed  that  sufficient  notice 
is  given  so  that  all  major  120 -hour  checks  will  be  completed  prior  to 
the  operation.      This  assumption  can  be  modified  with  an  appropriate 
adjustment  in  the  mean  repair  rate.      The  concept  of  maintenance  crews 
assigned  to  hangar  deck  areas  ("spots"),    as  developed  by  Collins  [3], 
will  be  used.     Each  crew  will  be  capable  of  all  types  of  maintenance 
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and  will  operate  independently  at  the  identical  mean  repair  rate    X. . 
The  number  of  spots  is  determined  by  the  average  number  of  such 
crews  available  to  work  continuously  around  the  clock  on  a  watch  basis. 

The  state  of  each  a/c  is  assumed  to  be  statistically  independent  of 
that  of  others,    and  the  launching  and  landing  transition  probabilities 
will  be  developed  on  the  basis  of  independent  Bernoulli  trials.      The 
parameters  can  be  determined  using  the  maximum  likelihood  estimators. 
The  range  of  the  number  of  a/c  desired  on  station  at  any  given  cycle  will 
be  set  by  the  user.     The  number  to  be  launched  at  any  time  is  assumed 
equally  likely  within  this  range.      This  input  parameter  is  a  function  of 
the  estimated  submarine  density  (i.  e.  ,    expected  contact  rate).      The 
lower  limit  will  be  set  at  the  number  of  a/c  desired  on  station  in  the 
search  (screening)  phase,    and  the  upper  limit  is  set  at  the  maximum 
practicable  number  of  a/c  to  be  launched  during  a  multiple -contact 
phase. 

Briefly,    the  assumptions  are: 

1.  a/c  will  be  relieved  on  station. 

2.  Any  desired  length  of  operation  can  be  set  as  an  input. 

3.  Major  120-hour  checks  will  be  completed  prior  to  the  operation. 

4.  No  resupply  to  the  carrier  is  available. 

5.  The  launch-to-launch  cycle  for  all  ASW  a/c  is  four  hours. 

6.  Minor  maintenance,    refueling,    and  rearming  only  can  be 
performed  on  the  flight  deck. 
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7.  Each  maintenance  spot  is  characterized  with  an  independent 
exponential  repair  time  with  mean  repair  rate  of   X.   for 
around-the-clock  operations. 

8.  The  number  of  a/c  lost  due  to  attrition  is  a  Poisson  random 
variable  for  each  cycle  period  with  parameter    A.       (a/c 
accident/flying  hours  for  a/c  type). 

9.  Any  a/c  lost  by  accident  will  not  be  returned  to  service  due  to 
either    (a)  physical  loss  at  sea,     or    (b)  insufficient  mainte- 
nance capability  aboard  ship  and  lack  of  major  parts. 

10.       The  number  of  a/c  launched  for  each  cycle  is  uniformly 

distributed  between  the  upper  and  lower  limits  determined  by 
the  user. 
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3.       MODEL  DESCRIPTION 

3.  1    The  Theory 

A  stochastic  or  random  process  is  a  collection  of  random  variables 
indexed  on  some  set  T,    (X(t),    te  T).      In  this  case,    time  is  the  indexing 
set,    and  the  Markovian  assumption  states  that  the  future  state  of  the 
process  depends  only  on  the  state  at  the  present  time  and  not  on  its 
past  history.      Due  to  the  cyclic  nature  of  our  problem,    it  is  possible 
to  increment  time    (T  =  (0,    1,    .  .  .  )  )  using  the  cycle  time  from  launch 
to  launch  as  the  steps  of  unit  time  in  a  discrete  Markov  chain.      It  is 
assumed  that  the  reader  is  familiar  with  the  notion  of  a  random  variable 
as  a  function  defined  on  a  sample  description  space  (S)  on  which  the 
family  of  events  or  outcomes  (E)  of  a  probability  function  can  be 
defined    [4], 

A  discrete  time  Markov  chain  is  described  by  a  sequence  of  dis- 
crete valued  random  variables  and  is  determined  when  the  one -step 
transition  probabilities  of  the  state  variables  are  specified,    i.  e.  ,    a 
conditional  transition  probability  of  a  transition  at  time    n   for  each  pair 
of   i,  j  =  0,   1 ,  .  .  .   ,  m   (m  being  the  number  of  states  in  the  process) 
must  be  given. 

p..(n,   n+l)  =  P[X(n+l)=j|X(n)  =  i] 
If  the  transition  probability  functions  depend  only  on  the  time  difference, 
we  have  time  homogeneity 
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p     (n+  1,1)    =    p..(0,   1)     =    p..    . 

The  initial  state  of  the  system  must  be  given  either  as  a  specific  state 
or  randomly  as  a  probability  distribution  function  over  the  possible 
states. 

The  p..  (transition  probabilities)  are  arranged  in  matrix  form  and 
satisfy: 

1.  p. .  s  0      for  i,j  =  0,   1 ,  ...   ,m; 

ij 

m 

2.  £      p..    =     1  ,    i.  e.  ,    the  rows  of  the  transition  matrix  sum  to  1 

for  all  i  for  the  states  within  the  description  space  [4]. 
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3.  2    The  Model 

In  order  to  establish  the  finite  set  of  states  (E)  for  the  model,    we 
shall  consider  two  random  variables  defined  as  follows: 

X    (t)     =     The  number  of  a/c  flying  at  time  t   not  having  flown  in 

the  previous  launch-to-launch  interval. 

X    (t)    =     The  number  of  a/c  in  or  awaiting  maintenance  at  timet. 
c* 

Now,    we  will  consider  the  vector    X(t)  =  [X    (t)  ,    X    (t)  ]  as  a  pair  of 

random  variables  and  thereby  have  a  bivariate  stochastic  process  with 

the  possible  states  ranging  from  (0,  0)  to  (A,  N). 

0  ^  X    (t)  ^  A    =    No.    of  a/c  desired  on  station,    and 

0  <•  X    (t)  ^  N    =    No.    of  a/c  of  given  type  aboard  carrier. 

We  will  define  an  operating  cycle  as  an  interval  unit  of  time.      Process 
observations  of  X(t)    will  be  made  at  successive  unit  interval  launch 
times.      To  develop  the  p..  elements,    consider  a  given  time  t   for  launching 
until  A  aircraft  are  flying  or  until  the  supply  of  ready  a/c  is  depleted. 
Those  a/c  failing  the  launch  enter  the  maintenance  state  at  this  idealized 
point  in  time  t  (the  total  launching  time  required  is  much  less  than  the 
total  cycle  time).     At  some  time  T,    less  than  the  launch-to-launch  unit 
time  interval,    the  a/c  which  were  relieved  on  station  return  and  land  at 
the  idealized  point  in  time   t  +  T.     Some  of  these  a/c  will  require  mainte- 
nance and  enter  the  maintenance  queue.      Those  requiring  only  refueling 
and  preflight  inspection  will  enter  a  ready  status  to  be  tested  for  the  next 
launch. 
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During  the  unit  time  interval,   maintenance  will  be  performed  on 
those  a/c  in  the  not-ready  status,    and  a  certain  number  of  aircraft  will 
be  repaired  according  to  assumption  7. 

In  summary,    we  start  the  system  in  some  initial  state  (such  as 

(0,  0)  with  no  a/c  flying  or  in  maintenance)  or  start  with  a  probability 

distribution   Q(t   )  over  the  states,  E,    at  time  t    .      We  launch,    recover, 
o  o 

and  repair  a/c  in  the  unit  interval  and  repeat  the  process  over  each 

succeeding  unit  time  interval  until  the  end  of  the  operating  period. 

Knowing  the  transition  probabilities  within  the  unit  time  interval,    we  can 

develop  the  elements  of  the  transition  matrix,    P,    or  {p  .       ..       ._     ..  }  . 
*  ltr(a,  i)  ,    (p,  j) 

These  are  the  probabilities  of  going  from  the  state  of   a   a/c  flying  and 
i   a/c  in  maintenance  to    {3    a/c  flying  and  j    a/c  in  maintenance  over  the 
unit  time  interval. 

It  was  assumed  in  section  2  that  A,    the  number  of  a/c  to  be  launched, 
and  N,    the  total  number  of  a/c  on  board,    are  random  variables,    whereas 
they  have  been  treated  as  constants  so  far  in  the  development.      To  be 
analytically  correct  in  including  this  feature,    one  should  develop  the 
appropriate  quadrivariate  process.     Such  a  development  leads  to  too 
large  a  state  space  and  the  author  chose  to  include  these  effects  by  using 
a  Monte  Carlo  simulation  technique.      That  is,    at  the  beginning  of  each 
cycle,    a  random  mechanism  is  used  to  determine  the  values  on  A  and  N. 

The  probability  of  losing  an  a/c  or  changing  the  desired  number  to 
be  launched  is  determined  from  the  specified  distributions  at  the  beginning 
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of  each  unit  interval,    and  the  resulting  P  matrix  containing  the 

p,  .  .    is  then  recomputed.      The  probability  distribution   Q(t) 

(a,  i)  ,    (P  ,  j) 

over  the  states  at  any  time   t   may  be  determined  by  the  appropriate 
number  of  successive  iterations  of  the  Q  vector  times  the  P  matrix,    i.  e.  , 
Q(t)    =    P[Xx(t)  =  0,  X2(t)  =  j]    =    Q(t  -  l)xP    . 
The  probability  of  maintaining   a   a/c  on  station  over  any  given  period 
of  operation  may  be  obtained  at  any  unit  time   t    (i.  e.  ,    the  beginning  of 
the  next  cycle)  by  summing  out  the  appropriate  maintenance  state  proba- 
bilities.     Thus,    P(a   a/c  are  flying  at  time   t  )  = 

Pr  (X    (t)  =  a  )     =  jj     Pr  (Xj  (t)  =  a  ,    X,,  (t)  -  i  )    . 

i=o 

A  mathematical  comment  appears  to  be  in  order.      In  the  case  of 

fixed  A  and  N,    the  states  of  the  Markov  chain  are  positive  recurrent;  and 

steady-state  probabilities  can  be  found  for  the  entire  state  space.      In  the 

case  of  decreasing  N  due  to  a/c  attrition,    this  is  not  true;  and  (0,  0) 

becomes  an  absorbing  state  as  time    (t)    goes  to  infinity.      This  latter 

consideration  is  not  a  realistic  one  for  the  operational  period  envisioned. 

Therefore,    it  is  mathematically  more  feasible  to  use  the  former  chain  in 

conjunction  with  the  Monte  Carlo  technique. 
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4.       DEVELOPMENT  OF  THE  TRANSITION  MATRIX 

Perhaps  the  simplest  way  to  view  this  development  is  to  note  the 
various  transition  probabilities  incorporated  in  one -unit  time  cycle 
defined  as  follows: 

(1)  Y  r  i     -   the  launching  transition  probabilities  at  time  t.      This  is  the 

fgh 

probability  of  taking  f  ready  a/c,    launching  g  successfully,    and  sending 

h  into  maintenance.      Each  a/c  to  be  launched  is  considered  a 

Bernoulli  trial  with  probability  of  failure  of  p    ,    which  is  estimable 

Y 

and  subject  to  sensitivity  analysis.      The  values  of  y         are: 

fgh 

a.  0      if  g  >  A  ,      since  only  A  a/c  are  desired; 

b.  0      if  g  +  h>  f;    it  is  impossible  to  launch  and  send  into  mainte- 
nance more  a/c  than  are  available; 

c.  0      if  g  <  A,    g  +  h  <  f  ;    launching  continues  until  A  a/c  are  flying 
or  until  all  f  are  used  up; 

f  g  f  -  a 

d-       („)  (!    -  P    )     (P     )  if  g  <  A,    g  +  h  =  f ,     standard  binomial 

g  v  Y 

when  all  a/c  in  the  ready  state  are  used  up  but  the  A  a/c  are 

not  launched; 

g+  h  -  1  g         h 

e«       (        h       )    (1   -  p)      (p)         if  g  =  A,    g  +  h  >  f ,    standard  negative 

binomial  for  g  successes  in  g  +  h  -  1  trials. 

(2)  II     (m)    =   the  landing  transition  probabilities  which  occur  at  time 
t+  T  .      We  must  consider  the  probability  that  if  there  are  a/c  flying 
at  time  t  then  m  a/c  will  enter  maintenance  at  recovery  time  t+  T. 
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II     (m)  will  equal  a  standard  binomial  where  p  =  the  probability  of 

a 

equipment  failure  in  flight: 

n    (m)   =    (a  )  (1   -  pf  "m(p)m  ,     m=  0,  1,  ...  ,a  . 

a  m 

(3)     P..(t)  =  the  maintenance  transition  probabilities,    i.  e.  ,    the  proba- 
bility of  repairing  (i  -  j)    a/c  in  time  t  .      Two  maintenance  periods 
occur:  the  first  starting  at  time  t  and  ending  at  time  t  +  T  ,    the 
second  starting  at  time  t  +  T  and  ending  at  the  end  of  the  cycle, 
(t  +  1).      Under  assumption  7,    the  pulsed  input,    multiple  exponential 
server  queue  is  developed  with  D  maintenance  "spots"  or  servers 
each  with  identical,    independent  service  rates,    X  .      For  each  server, 
then,    the  probability  of  remaining  occupied  (given  the  server  is 

-  X  T 

busy)  in  time  t  =  e  .      The  probability  of  becoming  free  (i.  e.  , 

repairing  an  a/c)  =  1    -  e  .      The  resulting  queueing  equations  are: 

A.  dP.        (t)  /  dt  =    -  n  \  P.       (t)   +   (n  +  1 )  X  P.  ,  (t)      for  0  <:  n  <  D  ; 

l,  n  i,  n  l,  n+  1 

B,  dP.        (t)  /dt    =    -DXP.        (t)   +    DXP.  .  (t)  for  n  s>  D   . 

i,  n  i,  n  1,11+1 

Three  ranges  of  i   (initial  queue  state)  ,     j    (final  queue  state)  ,    and 

D  become  significant: 

a.       When  j  <.  i  <.  D,    then  not  all  spots  are  busy  since  there  are 
fewer  a/c  in  maintenance  than  spots.      Each  spot  works  inde- 
pendently; therefore,    the  solution  to  A  is  the  binomial: 

(i  -  j) 


- \t.   x        "         -  Xtj 


P;;(t)       =         ())       (1      "     e^        ) 
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b.  When  D  ^  j  ^  i,    then  all  spots  are  occupied  throughout  the  total 
service  time,    and  the  closed  form  solution  to  B  is  the  Poisson: 

pij(t)  -  (TTiyi 

c.  When  j  <  D  <  i ,    then  all  spots  are  busy  at  the  beginning  of  the 
service  period,    and  some  spots  become  idle  during  the  service 
period.      The  explicit  form  solution  of  equation  A  is  found  using 
moment  generating  function  transformation: 

Pyw  -  V  <?><?>  {(o^r)  (i-D,e-xtn 

k=o 
(The  derivation  of  this  solution  is  discussed  in  Appendix  I.  ) 

The  figure  on  the  following  page  will  show  the  relationships  of  these 
transition  probabilities  within  the  unit  time  interval. 
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a/c    ON  STATION 


LAUNCH   a/c 


ffgh 


READY   a/c 


pi,j(T) 


LAND   a/c 


g-m 


il     (m) 
g 


m 


pk>i(1-T> 


MAINTENANCE 


m 


TRANSITION  PROBABILITIES  WITHIN  THE  UNIT  CYCLE 


FIGURE  1 
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In  order  to  develop  each  transition  probability  over  the  total  unit 
time  interval,    we  must  consider  all  events  taking  place  within  the 
interval.      Thus,    to  obtain  the  probability  of  going  from   a   a/c  flying  and 
i   a/c  in  maintenance  to    (3    a/c  flying  and  j   a/c  in  maintenance,    we  start 
at  the  state  (a,  i)  at  time   t.     At  this  time,    a/c  are  launched  and  some 
1    a/c   failing  the  launch  enter  maintenance.      These    i  +  1   in  maintenance 
are  then  serviced  until  time   t  +  T   when  some   k   a/c  are  still  in  the 
maintenance  state.     At  time   t  +  T,    of  the    a    a/c  previously  flying,    some 
m  enter  maintenance  and  (a  -  m)  enter  the  ready  pool.     Maintenance  is 
continued  on  the  (k  +  m)  a/c  for  the  remainder  of  the  cycle  (1   -  T)  ,    until 
the  end  of  the  unit  time  interval  when  j    a/c  remain  in  the  maintenance 
state.      In  functional  form: 

N-a-i      i+1        a 

P(«,  i),    (3,  j)    =      ^        VL         L        YN-a-i,    B,    1 

1  =o       k=o     m=o 


Pi+l,k(T>-na(m)    -Pk  +  m.j'1-7' 
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5.       SUMMARY 

Representative  values  for  the  mean  repair  rate  and  the  landing  and 
launching  failure  rates  produced  results  in  agreement  with  the  sensitivity- 
analysis  by  Collins  on  these  parameters  in  [5].     For  failure  proba- 
bilities less  than  .  5,    and  mean  repair  rate  less  than  12  hours,    the  effect 
of  reducing  the  available  maintenance  time  to  80%  of  the  cycle  time  was 
negligible.      Optimal  loading  and  cycling  policies  can  be  determined  for 
known  values  of  these  rates. 

The  model  affords  the  following  checks:    (1)  the  rows  of  each  P 
matrix  are  summed  as  they  are  computed  by  the  program;  and  (2)  the 

probability  distribution  vector  (QJ)  is  summed  over  the  states.     Each 

-  8 
summation  was  within  10  of  one  in  the  computer  model. 

The  user  may  substitute  any  available  distribution  over  the  interval 

of  a/c  desired  on  station.      In  order  to  keep  A  fixed,    enter  the  desired 

value  as  both  upper  and  lower  limit  (  A  =  ALOLIM  =  LUPLIM  )  .     For 

-  8 
fixed  N,    use  a  very  small  value  for  ALAM  (such  as  10         )  .     Subroutine 

KRAN  is  a  uniform  generator,    using  the  half  open  interval  (lower  limit  +  1, 

upper  limit  +  2)  and  a  starting  number  as  inputs.     KRAN  outputs  an  integer 

in  this  interval.     Subroutine  DRAW  was  used  to  provide  some  intuitive 

grasp  of  the  results.      DRAW  was  used  in  binary  card  form  and  is  not 

essential  to  the  main  program.     (The  indicated  associated  statements 

must  be  removed,    however.  ) 
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The  results  of  reasonable  arbitrary  parameter  values,    based  on  the 
author's  experience,    have  shown  that  most  of  the  probabilities  concen- 
trate over  a  few  states.     Moreover,,    computation  time  increases  rapidly 
as  a  function  of  N  (no.    of  a/c),    see  Figure  2.      This  would  indicate  that  a 
simple  approximation  to  the  model  could  be  developed.      One  method 
presently  being  investigated  to  reduce  computation  time  is  to  shrink  the 
probability  state  space  to  include  only  those  significant  states  and,    thus, 
reduce  the  size  of  the  transition  matrix.     Alternatively,    the  eigenvector, 
eigenvalue  representation  of  the  P  matrix,    might  be  used. 

Originally,    it  was  hoped  to  utilize  the  data  from  the  Fleet  ASW  Data 
Analysis  Program  ( FADAP)  to  attempt  a  verification  of  the  model  with 
its  real-world  counterpart.      The  only  method  available  at  present  for 
obtaining  the  necessary  data  is  by  direct  observation  or  a  program  of 
data  collection,    as  suggested  by  Collins  [5]  . 

Many  fruitful  areas  of  investigation  exist: 
(1)     Attrition  has  been  simply  modeled  by  the  Poisson  method.      The  two 
components  of  attrition,    accidents  and  supply  shortage,    can  be  more 
accurately  modeled  and  used  to  develop  logistic  schedules  for 
maintenance  and  supply.      One  simple  technique  is  to  assume  each 
component  is  independent  and  Poisson,    and  estimate  a  supply  failure 
rate  for  AOCP  attrition  from  past  data.      With  these  assumptions, 
the  total  attrition  is  Poisson,    with  the  parameter  equal  to  the  sum 
of  the  accident  and  supply  failure  rates. 
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(2)  The  model  can  be  modified  to  make  the  number  of  maintenance  spots 
available  for  any  cycle  a  variable  function  of  time,    D(t). 

(3)  An  investigation  of  the  Markovian  assumption  validity  as  the  cycle 
times  become  smaller  and  smaller. 

(4)  Development  of  a  continuous  time  model. 

(5)  Modification  of  the  model  to  simulate  re  supply  by  COD. 

(6)  A  study  of  the  distribution  of  submarine  contacts  to  determine  the 
validity  of  the  uniform  a/c  demand  assumption. 
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APPENDIX  I 

EXPLICIT  SOLUTIONS  OF  THE 
MAINTENANCE  QUEUEING  EQUATIONS 

The  queueing  equations  for  the  pulsed  input  queue  are  essentially 
the  pure  death  process  given  in  [  1  ]   and  [4]   as  problems  and  developed 
by  Collins  in   [5].      The  equations  are: 

dPi     n(t) 

A.      i- =     -n\P,        (t)    +    (n+l)\P.  ,  (t)     for  0  <  n  <  D 

dt  l,  n  l,  n  +   1 


dPi    n(t) 

B.      j- =     -nXP.        (t)    +     D\P.  ,  (t)  for  n  ^  D 

dt  l,  n  i,  n  +   1 


where    P. .  (0)  =  A. .    and  P. .  (t)  =  0  for  i  <  j  ,    since  no  input  (arrivals) 
ij  ij  ij 

occur  during  the  service  time. 

Equation  B  is  solved  directly  in  closed  form: 

(n  -  i)  :    -  \Dt 
P.        (t)     =    {DXt) 


i,  n  (n  -  i)  ! 

Now  transforming  the  first  equation  (A)  using  the  moment  gener 

ating  function  (MGF),  g 

D-l 
G(s,  t)    =      £      s    P    (t)    , 
n=0 

as  outlined  in  [  4  ]    (Chapter  7),    and  its  partial  derivatives: 
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,1)  .£-  =  D£'   snP'(t) 

dt  n 

n=0 


(2)  4^    =         ^       nsn-   1  p 

dS  n=0 


Where  P    (t)  denotes  the  conditional  probability  P.        (t)  ,    by  substituting 
n  i,  n 

(A)  into  (1),    properly  identifying  the  first  summation  with  (2),    and 
changing  the  second  summation  index  to    r  =  n  +   1  ,    we  get: 


dG  dG  ,2 


X  s   — —     +     X       Z      r  s  P    (t)    ,  or 


dt  ds  r 

r  =  0 


<3>  f^^'-'i't*  lD'D"lpDf"    • 


since 

D 


rs'-'P(t)  ,«  +  Di»-'pd„, 

r  =  0 


Next,    replace  the  partial  differential  equation  (3)  with  a  system  of 
ordinary  differential  equations  using  the  Lagrangian  auxiliary  equations: 

dt  ds  dz 


1  ^<S-1»  kDSD-'V„ 


The  solution  to  the  first  equation  (using  the  first  two  differentials) 
is: 

Xt    =    ln(s-l)    +     C' 
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and  hence 


or 


s     =     C    e  +     1 


Gx     =     e~Xt  (s   -  1)     . 


The  second  equation  is:    (using  first  and  third  differentials) 


dz    =     -  \D  (C    e  Xt   +    1)°  '  1  P     (t)  dt    . 


Using  the  solution  to  (B)    where  m  =  i  -  D  to  replace  P     (t)    and 
integrating,    term  wise,    the  binomial  expansion  of  (C    e        +1) 


z=UD)m+1    Dr    rD:1)   CJ    ft».-MD-J>tdi 


^crx^.i' 


where  the  integral  is  evaluated  as: 

m  k     -\  (D  -  j)t 

■o  t    e  m  „ 

k=o    (MD-j))m-k+1       k!  2 

Thus, 

*  t-w     D_1  m        /  ,  r^.vk  «    v    m  -  k 

C       =     z+e"XDt  'D^    -       -J      -        (XDt> 


j=0       J  k=0         k!         ^D-iJ 


and  the  general  solution  is   0  (C    ,    C    )  ,    where  0    is  an  arbitrary  function 
and 

C       =    u  (s,    t,    z) 

and 

C^     =     v  (s,    t,    z)     . 
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To  get  our  particular  solution,    use  the  boundary  conditions  for    G(s,  t) 

(1)     for  s  =   1  , 

D-l 

G(i,  t)  =    r    p  (t) 

n=0         n 
=    Pr  [no.    in  maintenance  at   t    is  <  D  I    i   at   t  =  0  ] 


i-D  "^t  n 

G(l,    t)    =    1     -       Z       ~ n,UDt)         =,■!«-♦    (t) 

n=0 

where 

u  (1,   t,    z)    =     C       =    0 

fi'V      I          r                 j    »-XDt     ^       (XDt)k 
v  (1,    t,    z)    =     C2    =    z  +  e  L     — 

k=0 
so 


C2    =    z  +  tx(t) 

(2)     for   t  =  0  , 

D-l 

G(s,  0)     =      £        snp    (0)     =    0    ,       since    i  ^  n>  D 

n 
n=0 

where 

u  (s,    0,    z)     =     (s   -  1) 

D-l  .  n         m 

v(s,    0,    z)     =     C2     =     z    +      Z        (  j)     (s   -  1)J    (5-7-) 

Thus, 


D-l  m 

G(s,    0)     =    z    +      £        (?)     CJ     (f^-^— )  -     C. 

j=0         J  VD-jy 
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Substituting  the  general  value  for  C      above  : 

D-l         _  ■         .  m 

G(s,    t)     =    0  (u,    v)     =      2        (J)    (s  -  1)J  e  _AtJ     C^~) 

3  =  0  VD-jy 


j=0  k=0  V        "  jy 


Rearranging  terms, 

D-l  D-l  _  .  m 


(s,  t)  =  l  sn  e  (jj  <? )  (-dj  r  r^^ 

n=0  j  =  n       n  J  L  V.  D  -  j^ 


-  \tj 
e 


-XDt      m 
-  e 


k=0 


where    P    (t)     =    the  coefficient  of   s 
n 


»         UDt)K      ^      D    Nm-K     1 
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APPENDIX  II 
THE  LOGICAL  FLOW  DIAGRAM  OF  THE  COMPUTER  PROGRAM 
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SET  ?W*.   TRANSITION  MATRIX  =0 
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COMPUTE  p,    w    At)   =  ?R(I,IBETA,IY) 

; \      ?XJ  \      iJ  J 1 

i 

(  I  =  IX+(lALPH-l)xTOTAL  A/C ) 


PRINT  THE  p^-'s  OF  THE  (IALPK,IX)TH  ROW  AND  ROW  SUM 
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i 
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OR 

=0 
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1 
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QI 

7 

X 

P 
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COMPUTE  THE  (IBETA,  IY)TK  ELEMENTS  OF  QJ(J) 


(J=IY+(IB.TA-l)xN) 


SUM  THE  ELEMENTS  OF  THE  QJ  VECTOR 


, _j 

[COMPUTE  THE  MEAN  OE  A/C  LAUHGHEdI 

■r~ — i 1 


l 
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PLOT   QJ(J)    VS    J 


-viaa. 


!SET  QI  =  QJ 


. 


T 


605  Hs- 


REINDEX  MATRIX 
llALAST  «IA 
INLAST  =  NEWN 

IALASTM1  =  IALAS' 


STOP 
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APPENDIX  III 
THE  COMPUTER  PROGRAM 
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:B?NAR^56?AN,0/49/S/1S/25/E/45  =  54,15  '30000'5.  000 

(RFLOCOM.  000 

-FTN.E.  00° 

PROGRAM  MARKOV  °0° 

OPPPA??n^AM  ILA    ^STATIONARY  BIVARIATE  MARKOV  CHAIN  MODEL  OF  ASW  A/C  000 

OPERATIONS.   THE  RANDOM  VARIABLES  ARE  THE  NUMBER  OF  A/C  FL YINfi  AT  Tnp  aaS 

thVZS?    %    ANY  GWEN  LAUNCH  CYCLE-   THE  MAX?mSm  NO*  OF  A/C  ALlJweS'in  000 

F    ; ELJS  16(NA>-   THE  RANGE  OF  A/C  TO  BE  LAUNCHED  AT  ANY  GIVEN  000: 

in  lrJ    In    6  A/C-   ™E  BLOWING  INPUTS  ARE  REQUIRED.  qoI 

ID=  THE  NO.  OF  INDEPENDENT  MAINTENANCE  SPOTS  on? 

NA=  TOTAL  NO.  OF  A/C  TYPE  ON  BOARD  ™ 

I.'am'mIam  ^2M  LAUNCH  T°  RFCOVERY/LAUNCH  TO  LAUNCH  CYCLE  TIME(HRS)001 
FLAM=MEAN  REPAIR  TIME  PER  SPOT/LAUNCH  TO  LAUNCH  CYCLE  TIME     HRsSo' 

P- ARSn   S^Jt1!:  °F  A/C  FATLING  LAUNCH(M.L.EST.  FROM  PAST  DATA)  00 1 

?M  A/Cp^t!LURE  DURING  FLIGHT  QUIRING  MAINTENANCE  00 1 

AT  LANDING  (M.L.  ESTIMATOR  FROM  PAST  DATA)  nn i 

01=  THE  PROBABILITY  DISTRIBUTION  VECTOR  OVER  ALL  POSSIBLF  STATES  00 1 

7  X  17  =  119)  SUCH  THAT  THE  SUM  OF  ALL  QI ( I )  =  l.    THIS  001, 

TCYMF    ToJE°r,rLlHl    cSER  AND  INPUTTED  BY  USING   A  DATA  STATEMENTS  V 

ICYCLF  =  NO.  CYCLES  DESIRED  FOR  OPERATION  no?r 

JCYCLE  =  LAUNCH  TO  LAUNCH  T IME ( HRS ) ( TOT.  TIME=ICYCLE  X  JCYCLFi  • 

ALAM  =  ACCIDENT  RATE  FOR  TYPE  A/C  (ACCIDENT/HOURS)  °0°0l 

ALOLIM  =  DESIRED  LOWER  LIMIT  ON  A  ~nl 

AUPLIM  =  DESTRED  UPPER  LIMIT  ON  A  °J' 

COMMON  FLAM,TIMF  002< 

TYPE  DOUBLE"  FLAM  002' 

COMMON  PTRM,GAMMA,PR,PRFMA,ID  °02' 

DIMENSION  BC(17),A(17),FBC(17)  002' 

DIMENSION  PTRM(17,17,2),GAMMA(17,7,17),PRFMA(7,7)  2?] 

DIMENSION  PR(119,7,17),QI(119),QJ(119)        '  °?1 

DIMENSION          FJPL0T(119),JT(12)  °2?f 

ENTER  DATA  CARDS  HERE  00ZS 

DATA((QI(I),I=1,119)=.2,16(.05)»102(.0))  003C 
NA  =  16 

ALAM=.01 

10  *  8 

FLAM=3.0 
PGAM=P=.4 
IYY  =  13421773 
TIME  =  .125 
ICYCLE=20 
JCYCLE=4 
AL0LIM=4. 
AUPLIM=6. 

END  OF    DATA  CARDS 

AL  =  AL0LIM+1.   $  AU  =  AIJPLIM  +2.  °°H 

UNITT=1.                                    •           '  0032 

N=NA+1  °033 

IAMAX=7  0034 

IALAST=0  °035 

D=FLOATF(ID)  0036 

NLASf=NEWN=N  0037 

KT=1                                                                                  N  90 3§ 

809     IA=KRAN(AL.AU,IYY)  °039 

IF(KT.EO.l)     113*115  0040 

0041 

42 


115  T1=-L0GF( .000000001  +  RANF (-1 ) ) *2. 30258/ALAM  0042 

IF(Tl-TFLC) 130,131,132  0043 

130  T2=-LOGF( .000000001  +  RANF ( -1 )) *2. 30258/ALAM  0044 
IF(T1+T2-TFLC)  230, 231, 1M  0045 

230  T3=-LOGF( .000000001  +  RANF ( -1 )) *2. 30258/ ALAM  0046 
IF(T1+T2+T3-TFLC)  331,331,231  0047 

331  NEWN=NLAST-3   $  GO  T0113  0048 

231  NEWN=NLAST-2   $  GO  T0113  0049 
132  NEWN=NLAST  $  GO  T0113  0050 

131  NEWN=NLAST-1  0051 
113  PRINT  8882,IA,NEWN  0052 

IF(NEWN-IA)  15,13,13  .  0053 

15  IA=NEWN  0054 

13  IF(IALAST)  11,12,11  0055 

12  CONTINUE  0056 
FROM  THIS  NEXT  STATEMENT  TO  NO.  483  IS  CONCERNED  ONLY  WITH  THE  GRAPH     0057 

DO  482  1=1,12  0058 

482  JT( I )=8H  0059 
JT(1 )=8HE(A/C)  =  0060 
JT(3)=8HSP0TS  =  0061 
JT(5)=8H  T  =  0062 
JT(7)=8HJ  VS  QJ  0063 
JT(8)=8HVECT0R  0064 
JT(9)=8H  N  =  0065 
JT(11)=8H  A  =  0066 
DO  483  1=1,119  <  0067 
FI=I  0068 

483  FJPL0T( I )=FI  0069 
IALAST=IA  0070 
DO  1235  1=1,17  0071 
DO  1235  J=1,IAMAX  "  0072 
DO  1235  K=l,17  0073 
GAMMA( I ,J,K)=0.0  0074 
THIS  PT  THE  LANDING  TRANSITION  PROBABILITIES  ARE  COMPUTED  0075 
DO  300  IAA=1,IAMAX  0076 
DO  301  MI=1,IAMAX  0077 
IF( IAA-MI )31»32»33  0078 

31  PRFMA( IAA,MI )=0.  0079 
.  GO  TO  301                                               N  0080 

32  MM1=MI-1  0081 
PRFMA( IAA,MI )=P**MM1  0082 
GO  TO  301  0083 

33  IAM1=IAA-1  . 0084 
"MM1=MI-1  0085 

BC(D  =  1.0  0086 

PROD=FLOATF( IAA-MI )  0087 

DO  50  IP=2,MI  0088 

AIP=FL0ATF( IP-1)  0089 

PROD  =  PROD  +1.0  0090 

50  BCUP)=PR0D*BC(  IP-D/AIP  0091 

IG0=IAA-MI  0092 

PRFMAI IAA,MI )=  ( BC ( M I ) * ( 1 .0-P ) ** ( IGO ) ) *P**MM1  0093 

361  eQNTfNUF  8§§* 

300  CONTINUE  0095 

AT  THIS  PT  THE  MAINTENANCE  TRANSITION  PROBABILITIES  ARE  COMPUTED  0096 

DO  100  IT=1,2  0097 
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IF(IT-1)25, 25,26 

25  TAU  =  TIME 
GO  TO  28 

26  TAU  =  UNITT-TIMF 
28  DO101  1=1, M 

DO  102  IJ=1,N 
IF  U-IJ)  14,199,17 
199  IF( I-ID)  19,19,1999 

1999  PTRM(I,Ij,lT)=EXPF(-FLAM*TAU*D) 
GO  TO  102 

14  PTRM(I,IJ,IT)=0. 

GO  TO  102 
19  FJM1=FL0ATF( IJ-1) 

PTRM(I,IJ,IT)=EXPF(-FLAM*TAU*FJM1) 
GO  TO  102 

17  IF(I-ID-l)  1,1,2 
1  BC( 1  )  =  1.0 

PROD=FLOATF( I-IJ) 

DO  10  IP  =2,IJ 

AIP=FLOATF(IP-l) 

PROD  =  PROD  +  1.0 
10  BC(IP)  =PROD*BC(IP-l)/AIP 

ELT=EXPF(-FLAM*TAU) 

PTRM  (I.U.IT)=BC(IJ)*(1.-ELT)**(I-IJ)*ELT**(IJ-1) 
GO  TO  102 

2  IF(IJ-l-ID)  22,24,24 
22    CONTINUE 

CALL  PID( I,IJ,IT) 

GO  TO  102 
24  D=FLOATF( ID) 

FLDT=FXPF(-D*FLAM*TAU) 
FACT  =  1.0 

A(  1  )=1.0 

MM  =  I-IJ 

DO  20  M=2,MM 

FACT=FACT+1.0 

A(M)=A(M-1)*FACT 

PTRM ( I »IJ,IT)=(D*FLAM*TAU)**( I-IJ ) *ELDT/A( I-IJ ) 

CONTINUE 

CONTINUE 

GO  TO  120 

CONTINUE 

IF(  IA-IALAST)  120,121,120 

IF(NEWN-NLAST)111,117,111 

IALPH=IALASTM1  $  GO  TO  801 

CONTINUE 

dS1^0!1^-!^  LAUNCHING  TRA^ITION  PROBABILITIES  ARE  COMPUTED 

IGM  =  XMINOF  (  IA,IFF) 
DO  203  IG=1,IGM 
IGM1=IG-1 
D9  262  IH=1»N 
IHM1=IH-1 

BPROD=(  (l.-PGAM)**IGMl  )  *(  PGAM*-*  IHM1  ) 
86  IF(IG-IA)  91,87,84 
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20 
201 
102 
101 
100 

11 

121 

117 

120 

AT 


OC 
OC 
01 
01 
01 
01 
01 
01 

oi: 
oi: 
oi: 

01 
01, 
01 
01. 

oi; 
oi: 
or: 
oi: 

013 

01] 
01] 
015 
012 
012 
012 
012 
012 
012 
012 
012. 

012; 

013 

013: 

013- 

013: 

013: 

013: 

013t 
013: 
013G 
013< 
014(. 

014: 
014; 

014' 
014< 

oi  m 

014* 
014" 
0146 
0145 
0156 
0151 
0152 
0153 


91  IF( IG+IHM1-IFF)  84,82*84  0154 

87  IF( IG+IHM1-IFF)85,85«84  0155 

84  GAMMA(  lFF»IGt  IH)=0.  0156 
GO  TO  202  0157 

82  BC(1)=1.0  0158 

PROD=FLOATF( IFF-IG)  0159 

DO  30  IP=2.IG  0160 

AIP=FLOATF( IP-1 )  0161 

PROD  =  PROD  +  l.O  0162 

30  BC(IP)=PROD  *  BC(IP-1)/AIP  0163 

IHM1=IH-1  0164 

TFMP=  PGAM**IHM1      •  0165 

TEMP1={ l.-PGAM)**IGMl  0166 

BPROD  =  TEMP*TEMP1  0167 

GAMMA { IFF,IG»IH)=BC( IG)*BPR0D  0168 

GO  TO  202  0169 

85  FBC(1)=1.0  0170 
PR0D  =  FL0ATF( IGM1-1  )  0171 
DO  40  IP=2, IH  0172 
AIP  =  FLOATF(IP-l)  0173 
PROD  =  PROD  +1.0  0174 

40  FBC(  IP)=PROD*FBC{  IP-D/AIP  0175 

GAMMA ( IFF,IG»IH)=FBC( IH)*BPR0D  0176 

202  CONTINUE  0177 

203  CONTINUE  0178 

204  CONTINUE  0179 

REMOVE  CARDS  FROM  HERE  TO  NO  999  IF  PRINT  OUT  NOT  DESIRED  0180 

PRINT  9»( ( ( I ♦IJ»IT»PTRM( I »IJ»IT) »IT=1»2) »IJ=1»N) *I=1,N)  0181 

9  FORMAT  (1H1/(2(6H  PTRM ( I  2 » IH , I  2 » IH , I  2 . 3H )  =  E14.5)))  0182 

PRINT  99 »( (  ( IFF,IG»IH,GAMMA( IFF.IG.IH) »IFF  =  1,N) » I G=l » I A ) » IH=1 »N )  0183 

99  FORMAT( 1H1/(2(7H  GAMMA ( I  2 » IH , I  2 » IH , I 2» 3H )  =  E14.5)))  0184 

PRINT  999»( ( IAA»MI ,PRFMA( IAA,MI ) »IAA=1»IAMAX) ,MI«1,IAMAX)  0185 

999  F0RMAT(1H1/(2(7H  PRFMA ( I  2 ♦ IH, I  2 »3H )  =  E14.5)))                         0186 

NOW  THE  TRANSITION  MATRIX  MUST  BE  ZEROED  0187 

111  CONTINUE  0188 

DO  899  J=l»119  0189 

DO   899  K=l»7  0190 

DO  899  L=l,17  0191 

899  PR(J»K,L)=0.0  0192 

START  COMPUTING  THE  ELEMENTS  OF  EACH  ROW*  I  =  I X+  (ALPHA  -  1)  X  TOTAL  A./C  0193 

DO  1000  IALPH=1,  IALAST  0194 

801  CONTINUE  0195 

DO  1100  IX=1»NLAST  0196 

COMPUTE  THE  P  ELEMENTS  OF  THE  IAPH»IX  ROW  AND  SUM  THE  ROW                 0197 

TSUM=0.  0198 

I=IX+( IALPH-1)*N  0199 

DO  800  IBETA=1»IA  >                      0200 

RSUM=0.0  0201 

DO  900  IY=1,NEWN  0202 

PR( I,IBETA,IY)=0.  0203 

ILIM=NEWN-IALPH-IX+2  0204 

PSUM=0.0  0205 

$UM*0«0  020^ 

SUML=0.0  0207 

DO  500  IL*1»ILIM  0208 

KLIM=IX+IL-1  0209 
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PRINT  888  ,     TSUM,IALPH,IX 


CONTINUE 
QJ( J)=QP 


021 
021 
021 


IXPIL  =  IX  +IL  -  1 
SUMM=0. 

DO  600  MI=1,IALPH 
SUMK=0. 

-   DO  700  IK=1,KLIM  ~f 

IKPMI  =  IK  +MI-1  "2I 

IF( IXPIL-NEWN)  701,701,700  «?i 

701  IF(IKPMI-NEWN)  702,702,700  «?r 

702  GAMH  =  GAMMA(  I  L  IM , I  BET A  ,  I L)  "ol 
PTRMH1  =  PTRM( IXPIL, IK,1) 
PRFMAH  =  PRFMA( IALPH.MI ) 
PTRMH2  =  PTRM( IKPMI, IV, 2) 
SUM  =  GAMH  *  PTRMH1  *  PRFMAH  *  PTRMH2 
SUMK=SUMK+SUM 
PSUM=PSUM+SUM 

700  CONTINUE  °pl 

SUMM  =  SUMM  +  SUMK  li' 

600  CONTINUE  "Jl 

SUML  =  SUML  +  SUMM 

PSUM2  =SUML  .  » 

500  CONTINUE 

RSUM=RSUM+PSUM 

PR(  I,IBETA,IY)=PSUM 
900  CONTINUE  c\%tX 

TSUM=TSUM+RSUM  rXvL 

800  CONTINUE  r,„ 

0236 


02  li 

02  It 
022( 
022: 
022; 
022: 


0221 
02  26 
022S 
023C 
0231 


0241 
0242 
0243: 
02  44 
0245' 
0246; 
0247- 
0248: 
0249 


888  FORMAT  (  7H  TSUM  =  E15.8.2I5)  0237 
1100  CONTINUE  .  ";j 
1000  CONTINUE                                                                ^39 

C      REMOVE   CARD  FROM  HERE  TO  889  IF  P  MATRIX  PRINT  OUT  NOT  DESIRED       0240 
DO  889  J=l,17 
DO  889  K=l,7 
D0889  L=l,17 
I=J+(K-1)*N 

889  PRINT  890,(PR(  I,LP,L)  ,LP  =  1,IAMAX) ,K,J,L 

890  FORMAT(7E14.5,2HJ=I2,5HK=l,A,?HL=I2) 
DO  898  1=1,119 

898  QJ( I )=0,0 
C  NOW   MULTIPLY   QI  AND  P  TO  GET  QJ 
805  PRINT  807,KT,IALAST, IA 
807  F0RMAT(1H1,13HQ  VECTOR  CASE  13///  15,15)  OPS! 

DO  802  IBETA=1,7  r\%\r 

DO  902  IY=1,17  Q253 

CAT  THIS  POINT  CALCULATE  THE  (IBETA,IY)TH  ELEMENT  OF  THE  QJ  VECTOR  0254 

J=IY+( IBETA-1)*N 
QP1=0. 
QP  =  0. 

DO  2001  IALPH=1,7 
DO  2201  IX=1,17 
I  =  IX+('IALPH-1  )*N 
0P1=QI  (  I )*PR(  It  IBETAtlY) 

QP=QP+QP! 

CONTINUE  21 zf 

0263 


0255 
0256 
0257 
0258 
0259 
0260 
0261 


0264 
0265 


46" 


PRINT  8882, IBETA, IY, J,QP  0266 

882  FORMAT(2I4,4H  QJ(I3,3H  )=  E14.8)                                         0267 

90?  CONTINUE  0268 

802  CONTINUE  0269 
CHECK  THE  SUM  OF  THE  Q  VECTOR  ■    0270 

QSUM=0.  02  71 

DO  808  J=l,119  0272 

808  0SUM=QJ( J)+QSUM  0273 

PRINT  8883, OSUM  0274 

3883  FORMAK6H  QSUM=   E15.9)  0275 

DO  333  I  ■  18,119  0276 

K  =  CI-D/17         •  0277 

FK=FLOATF(K)  0278 

FMEAN=  FK*QJ( I )+FMEAN  0279 

333  CONTINUE  0280 

TFLC=FMEAN*FLOATF( JCYCLE)  0281 

PRINT  335,FMEAN  0282 

335  FORMATf  17HMEAN  A/C  FLYING  =  F10.4)  0283 
STATEMENTS  FROM  THIS  POINT  TO  THE  CALL  DRAW  STATEMENT  REFER  TO   GRAPH   0284 

JT(2)=ICODE(FMEAN)  0285 

JT(4)=IC0DE(D)  0286 

FKT=FLOATF(KT)  0287 

JT(6)=ICODE(FKT)  0288 

FN=FLOATF(NEWN-l)  0289 

JT( 10)=ICODE(FN)  0290 

FIAA=FLOATF( IA-1)  0291 

JT( 12)=ICODE(FIAA)  .                                    0292 

CALL  DRAW( 119,FJPL0T,QJ,0,0,4H     , JT , 0 , 0 ♦ 0 , 0 ,0 ,0 , 8 , 8 ,0 ,L AST   )  0293 

FMEAN  =  0.  0294 

NEXT  WE  MUST  MULTIPLY  QJ  AND  P  TO  GET  QK  AND  SO  ON, . . ( QK+. . .N )             0295 

KT=KT+1  0296 

IF(KT-ICYCLE)  803,803,806  0297 

803  DO  804  1=1,119  0298 

804  01 ( I )=OJ( I )  0299 
IALASTM1=IALAST  0300 
IALAST=IA  0301 
NLAST=NEWN  0302 
GO  TO   809  0303 

806  STOP  06  0304 

END  0305 

SURROUTINF  PID(I,J,IT)  03  06 

COMMON  FLAM, TIME  0307 

COMMON  PTRM,GAMMA,PR,PRFMA,ID  0308 

TYPE  DOUBLE  BCBDCPROD  ,DID3  ,0  I  D4,DID5  ,DEXP  0309 
TYPE  DOUBLE  DAN ,DI Dl ,D I D2 ,SUM,DN ,ANM1 ,FAC ,COF, PSUM, PTR, FLAM, TAU ,D     0310 

DIMENSION  PTRM( 17,17,2)  ,BC(11) ,BDC(11)  0311 

DIMENSION  GAMMA( 17,7,17) ♦  PRFMA(7,7),  PR(119»7,17)                     0312 

D=FLOATF(ID)  0313 

IDP1=ID+1  0314 

IF( IT-I )25,25,26  0315 

25  TAU  =  TIME  0316 
GO  TO  2  8  0317 

26  TAU=  1 .-TIME  0318 
28  CONTINUF  0319 

IMDP1=I-ID  0320 

PTR=0.0  0321 
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DE 


A  TIME  AND  D  TAKEN  N  AT  A  TIME 


10 


20 


11 


201 


200 
103 

102 

101 


c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


PSUM-0. 

DO  200  NJ  =  J»  ID 
VELOP  N  TAKEN  J  AT 
BC(1 )=1.0 
PROD=FLOATF(NJ-J) 
DO  10  TP=2,J 
ATP=FLOATF( IP-1 ) 
PROD  =  PROD+1.0 
BC(IP)=PROD*  BC(IP-1)/AIP 
BDC( 1)=1.0 

PROD=FLOATF( IDP1-NJ) 
DO  2  0  IQ=2,NJ 
AIO=FLOATF( IQ-1) 
PROD=PROD+1.0 

BDC( IQ)=PROD*BDC( IQ-1) /A  I Q 
COF=BC( J)*BDC(NJ)*(-1)**(NJ-J) 
ANM1=FL0ATF(NJ-1) 

DAN=D/(D-ANM1 ) 
DID4=DFXP(-FLAM*TAU*ANM1) 
DID1=(DAN**( I- 1  DPI ) )*DID4 
SUM=0. 
DN  =  0. 

DO  201  K=1,IMDP1 
FAC=1. 
KM1=K-1 
PROD=0. 

DO  11  IK=1,KM1 
PROD=PROD+l. 
FAC=FAC*PROD 
IMIDK=I-ID-K 

SUM=( (FLAM*D*TAU)**KM1 )*DAN**IMIDK 
DN=DN+SUM 

DID3=DEXP(-FLAM*D#TAU) 
DID2=DN*DID3 
DID5=DID1-DID2 
PSUM=COF*DID5 
PTR  =PTR  +PSUM 
CONTINUE 

PTRM( I #J» IT)=PTR 
CONTINUE 
CONTINUE 
END 
FUNCTION  KRAN(A,B,IY) 


THIS  ROUTINE  RETURNS  AN  UNIFORMLY  DISTRIBUTED  RANDOM  INTEGER 


THIS  ROUTINE  RETURNS  A  INTEGER  RANDOM  NUMBER  .GE.  TO  A 

•  LT.  B 

A  =  BOTTOM  LIMIT  (INCLUDED)  FOR  THE  RANDOM  NUMBER 

B  =  TOP  LIMIT  (NOT  INCLUDED)  FOR  THE  RANDOM  NUMBER 

SET  IY  ONLY  ONCE  IN  MAIN  PROGRAM  FOR  EACH  SET  OF  RANDOM 

SOMg  GOOD  STARTING  VALUES  FOR  IY  FOLLOW 


/   FAC 


NUMBERS 


13421773 
33554433 
8426219 
42758321 


03 

03 

03 

03 

03 

03 

032 

032 

033 

033 

033 

033 

033, 

03: 

033 

03: 

03; 

03: 

03 
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03! 
03! 
03! 
03! 
03! 
03! 
03! 
03! 
03< 
03i 

o: 

0 

o: 
o: 

0: 

03 

o: 
o: 

o: 

0: 

0 
0 


0 
0 

o 


56237485 
62104023 
ANY  OF  THESE  MAY  BE  USED 

THIS  ROUTINE  MAY  BE  USED  IN  FORTRAN  60  OR  63 

IY  =  3125  *  IY 

IY  =  IY  -  ( IY/67108864)  *  67108864 
FY  =  IY 
-   KRAN  =  FY/67108864,  *  (B-A)  +  A 
RETURN 
END 

FINIS 
XECUTER. 


0378 
0379 
0380 
0381 
0382 
0383 
0384 
0385 
0386 
0387 
0388 
0389 
0390 
0391 
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VECTOR  CASE 

5 

5 

1 

1 

QJ 

(   1 

1 

2 

QJ 

(   2 

1 

3 

QJ 

(   3 

1 

4 

QJ 

(   4 

1 

5 

QJ 

(   5 

1 

6 

QJ 

(   6 

1 

7 

QJ 

(   7 

1 

8 

QJ 

(   8 

1 

9 

QJ 

I   9 

1 

10 

QJ 

C  10 

1 

11 

QJ 

t  11 

1 

12 

QJ 

1  12 

1 

13 

QJ 

I    13 

1 

14 

QJ 

:  14 

1 

15 

QJ 

:  is 

1 

16 

QJ 

:  16 

1 

17 

QJ 

'  17 

2 

1 

QJI 

'  18 

2 

2 

QJI 

19 

2 

3 

QJI 

20 

2 

4 

QJI 

21 

2 

5 

QJI 

22 

2 

6 

QJI 

23 

2 

7 

QJI 

24 

2 

8 

QJI 

25 

2 

9 

QJI 

26 

2 

10 

QJI 

27 

2 

11 

QJ< 

28 

2 

12 

QJI 

29 

2 

13 

QJ( 

30  : 

2 

14 

QJ( 

3i  : 

2 

15 

OJ( 

32  ) 

2 

16 

QJ( 

33  ) 

2  ' 

17 

QJ( 

34  ) 

3 

1 

QJ( 

35  ) 

3 

2 

QJ( 

36  ) 

3 

3 

QJ( 

37  ) 

3 

4 

QJ( 

38  ) 

3 
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APPENDIX  IV 
SAMPLE  RESULTS 

The  following  pages  present  the  values  of  the  elements  of  the  proba- 
bility distribution  vector  (Q  J)  and  its  graphical  plot  for  five  consecutive 
iterations,    i.  e.  ,    Q  x  P       for  n  =  1,    2,  .  .  .  ,  5  .      The  inputs  are  those 
shown  on  the  first  page  of  Appendix  III  between  statement  No.    30  and 
No.    31.      The  printouts  of  the  transition  matrices  and  their  computational 
elements  are  omitted.      The  plot  was  made  using  the  DRAW  subroutine  in 
the  U.  S.    Naval  Postgraduate  School  computer  facility  library.     Each 
vector  printout  contains  the  values  of  all  119  states  possible  (7  x  17)  and 
is  headed  by  the  past  value  of  A  +   1  and  the  next  value  of  A  +  1.     The  two 
indices  preceding  each  element  represent    3+1  and  j  +  1 ,    in  the  notation 
of  section  3.      For  example,    in  the  first  row  on  the  next  page,    the    "1     1" 
indicates  that  the  probability  of  being  in  state  (0,  0)  after  one  iteration 
is  =    .  026,    where  the  value  of  A  is  4  over  the  first  iteration.     Each  graph 
is  labeled  with  the  expected  value  of  a/c  flying,    the  number  of  mainte- 
nance spots  available,    the  vector  number  (T),    total  number  of  a/c  avail- 
able (N),    and  the  desired  number  of  a/c  on  station  (A).      The  "E"  notation 
indicates  the  power  of  10  to  multiply  by.      This  sample  run  demonstrates 
the  loss  in  total  a/c  and  variable  a/c  on  station. 
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